Different routes towards oscillatory zoning in the growth of solid solutions 
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Oscillatory zoning, i.e. self-formation of spatial quasi-periodic oscillations in the composition 
of solid growing from aqueous solution, is analyzed theoretically. Keeping in mind systems like 
(Ba,Sr)SC>4 we propose a ID model that takes into account the nonideality of the solid solution and 
the system asymmetry, in particular, reflecting itself in different solubilities for such systems. Based 
on a linear stability analysis different parameter regions can be identified. Even an ideal solution 
solution with a sufficiently large asymmetry can display oscillatory zoning. Numerical simulations 
complement the linear stability analysis as well as the qualitative consideration of the instability 
development and reveal the nature of the limit cycles. 

PACS numbers: 81.10.AJ, 47.54.-r, 05.65.+b, 82.40.Ck 
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I. INTRODUCTION 

Oscillatory zoning (OZ), i.e. spatial pattern made of 
quasiperiodic variations of the solid composition from 
core of crystals to their rim is widely met in natural 
minerals (see, e.g., Ref. The appearance of such 

patterns was traditionally related to cyclic changes in 
surroundings during crystal formation in rocks. How- 
ever the success of reproducing OZ in calcite crystals 0] 
and (Ba,Sr)S04 solid solutions [!, H], H[ in laboratory un- 
der quasistationary conditions has demonstrated the fact 
that at least partly OZ can result from selforganization 
during crystal growth in solution. 

The experimental setup used by Putnis et al. 0, 0, Q 
is sketched in Fig. [TJ It consists of two reservoirs, one 
filled with aqueous solution of BaCLo/SrC^ and the other 
with Na2S04. The two reservoirs are connected by a col- 
umn filled with silica gel to inhibit convective transport. 
At the beginning of experiments the reactants start to 
diffuse toward each other through the column. As the 
diffusion fields of Ba 2+ , Sr 2+ , and S04 2 ~ overlap and 
the solute concentration product exceeds the nucleation 
threshold in the vicinity of the column center, the crys- 
tal nuclei form. In approximately one month the experi- 
ments were terminated. The obtained crystals exhibited 
OZ although no external fluctuations were imposed on 
the system. 

Following the spirit of the general model by Ortoleva 
0, 0] for the growth instability caused by autocatalytic 
interaction of the species at the crystal surface L'Heureux 
et al. [8|, [9|, [lCj proposed a rather sophisticated model for 
OZ in the (Ba,Sr)S04 solid growing from aqueous solu- 
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FIG. 1: Experimental setup in which oscillatory zoned crys- 
tals of (Ba,Sr)S04 were synthesized by Putnis et al. 
The reactants counterdiffuse in the column and (Ba,Sr)S04 
crystals nucleate. The upper window sketches the structure 
of the nucleation zone and the length scales involved. 



tion. The detailed analysis of these models was carried 
out within the boundary layer approximation. 

In a previous paper [ll[ we have demonstrated that 
OZ in crystals growing from solution can be described as 
a boundary-reaction-diffusion problem. It is character- 
ized by passive diffusion of species through the solution 
bulk to the crystal surface where their interaction gives 
rise to the crystal growth. The latter, however, proceeds 
with a very low rate so that the crystal boundary can 
be treated as a surface fixed in space. It that work we 
have mainly studied the presence of the instability with 
respect to the nonideality parameter 9. It turned out 
that for sufficiently large 9, i.e. 9 > f? c o the instability 
and thus OZ can indeed be observed. 

Experimentally, however, it is observed that in partic- 
ular solid solutions which very different solubility prod- 
ucts of the endmembers display OZ (such as (Ba,Sr)S04) 
whereas systems with similar solubility products (such 
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as (Ba,Sr)CC>3) do not display OZ pjj]. For example, 
for the first case the solubility product of both endmem- 
bers differs by three orders of magnitude (see, e.g., [IH). 
The solubility product is related to the system asymme- 
try Thus, the question emerges whether the model 
also allows OZ for systems with a pronounced asymme- 
try rather than a significant nonideality. 

The purpose of this work is fourfold. First, we rederive 
our model in a somewhat extended way which allows one 
to better understand the microscopic origin of the differ- 
ent parameters. The definitions have been chosen such 
that the final model equations are identical to the model 
studied in our previous work [Tlj . Second, after deriv- 
ing the somewhat complex instability conditions from the 
linear stability analysis we argue on a semi-quantitative 
level that indeed the model possesses an additional in- 
stability channel for sufficient large values of the system 
asymmetry and thus obtain a semi-quantitative phase 
diagram of the instability region. Third, via a careful 
mathematical analysis we somewhat modify this picture, 
yielding some surprising features in the newly analyzed 
instability regime. Fourth, via numerical simulations we 
illustrate the behavior beyond the linear regime. 



II. MODEL 

A. Energetics of crystal growth 

We take into account the following mechanism of crys- 
tal growth having in mind the (Ba,Sr)S0 4 . The ions 
S0 2 ~ (below species of type 0), Ba 2+ (species 1), and 
Sr 2+ (species 2) diffuse to the crystal surface through 
the aqueous solution, where they are adsorbed and dis- 
play surface diffusion. If they reach the atomic steps 
they are incorporated into the crystalline lattice via the 
following precipitation reactions 

Ba 2+ + S0 4 2 ~ -> BaS0 4 (channel 0-1) , (la) 
Sr 2+ + S0 4 2 ~ -> SrS0 4 (channel 0-2) (lb) 

The latter process is considered to be irreversible, i.e. the 
solid dissolution is ignored, which means the system to 
be far from thermal equilibrium and the growth rate can- 
not take too low values in the case under consideration. 
Finally they are forming a new layer of the crystal. 

Migrating along the crystal surface adatoms experience 
many different local environments depending on the sur- 
face composition which will be characterized by the mole 
fraction \ of species 1 (0 < \ !)• The competition 
between adsorption and desorption is determined by the 
effective adsorption energies Ei(x) {i = 0, 1, 2) reflect- 
ing the species interaction with the crystal surface and 
aqueous solvent. In the mean field approximation they 



are written as 

Eq(x) = £o - 9iX -52(1 - X) , (2a) 
E 1 (x) = e 1 -g 1 +9(l-x), (2b) 
E 2 (x) = t2-g2 + 8 X , (2c) 

where all the energy quantities are measured in units 
of temperature, q is the solvation energy of species i, 
the constant gi characterizes the interaction between 
adatoms of type % = 1, 2 with atoms of type lying 
in the surface atomic layer of the crystal lattice, and the 
parameter 9 > quantifies the solid solution nonideality. 
This expresses the fact that the strongest interaction on 
the crystal surface holds between like ions. 

In these terms equilibrium between the adsorbed layer 
and the aqueous solution region adjacent to the crystal 
surface implies the following relation between the adatom 
concentrations Cj and the concentration Cf of the corre- 
sponding species near the crystal surface 

d = aCte- E ^ , (3) 

where a is the characteristic size of the crystalline cell. 

In principle, on vicinal crystal surfaces the adatoms 
should have some solvent shells and for them to be in- 
corporated into the crystal lattice these shells have to 
be destroyed. If it is essential then the precipitation re- 
actions (H|) at the surface atomic steps limit the crystal 
growth and the adsorbed layer can be assumed to be in 
quasiequilibrium, meaning equalities ([3]) to hold. In this 
case the partial rates $1 and #2 of the crystal growth 
though channels (jlap and (jlbp , respectively, are given by 
the expression 

a 6 

tfi = ^yc c,, (4) 

where Vi is the rate at which the pair of the Ba 2+ , S0 2 ~ 
adatoms or the Sr 2+ , S0 2 ~ adatoms meeting at the sur- 
face steps are incorporated in the crystal lattice and I is 
the mean distance between these steps. 

Combining expressions (O and ^ we get the desired 
relationship between the partial growth rates via the 
channels 0-1, 0-2 and the corresponding values of the 
solute concentrations Cq, Cf, and Cf near the crystal 
surface 

/,. \i/2 

i?i = w I -± J e-l" e^-e(i-x) C^Cf , (5a) 
/ \ 1/2 

$ 2 = w (ll\ e h e -*(i-x)-«x c s C s 2 . (5b) 
Here we have introduced the kinetic coefficient 

w = ^/I^ye 2 ^^ - 612 (6) 
and rewritten the interaction constants 51,2, £1,2 using 
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combination of the quantities 



1 



512 = - (.gi +92) , 
1 , 

£12 = - (ci + e 2 ) , 



= 5i - 92 , (7) 
?/ = ei - e 2 (8) 



to mark out the difference of species 1 and 2 in properties. 

Expressions (O are actually the main result of this sub- 
section and form the basis of the model for the crystal 
growth to be constructed in the next section. It is rather 
similar to the model we have developed previously [ll[, 
enabling us to sketch out the principle aspects only. Be- 
low we will assume the inequality > to hold be- 
forehand because, otherwise, the indices could be just 
exchanged. 



B. Model equations 

In the aqueous solution the SO4 - ions are assumed to 
be abundant. Thus, we can regard their concentration 
as a fixed value Cq. In this case the crystal growth in 
the ID description is governed by the boundary-reaction- 
diffusion model developed in our previous work [ll[. 
Namely, diffusion of the components i = 1, 2 through the 
solution is considered within the region z G [0, L] and is 
described by the equation 



dd(z,t) _ d 2 Ci(z,t) 



dt 



dz 2 



(9) 



where Di is the diffusivity of the species i in the aqueous 
solution and the system size L should be chosen large 
enough in order to enable us to fix the influx of both the 
components at the external boundary z = L 



Gi = D, 



dCi(z,t) 



dz 



(10) 



z=L 



Then having in mind expressions (O we write the follow- 
ing boundary condition at the crystal surface [z = 0) 



D, 



dd(z,t) 



dz 



z=0 



n(x) 



(ii) 



which relates the boundary values of diffusion flux and 
the rates of species attachment to the crystal surface 



aCf 



(12) 



caused by the growth process. Here the time scales of 
the crystal growth dynamics via the channels 0-1 and 
0-2 individually are specified as 



n(x) 



T 2 (X) = Tg 



1/2 

Yl\ P iv-<t>x+e{i-x) 



1/2 



.i r)+ 0(i_ x )+e x 



(13a) 
(13b) 



where the time scale of the crystal growth dynamics as a 
whole process is 



(14) 



Finally, the solid composition is governed by the equation 



dx 
dt 



(i-x) 



aC 9 s 



aCl 

n(x) x Mx) 



(15) 



following from mass conservation and used previously in 
a number of papers on OZ, see, e.g., Refs. H BBS El- 
The given system admits only one steady state solution 



G 

d(z) = Cf t +Xst-fT z > 

C2(z)=Cf, st + (l-Xst)^-2 



(16a) 



where G — G\ +G2 is the total diffusion flux determining 
the growth rate of the crystal as a whole, the correspond- 
ing value of the crystal composition x st = G\/G, so 



G\ — XstG , 



G 2 = (l-Xst)G, 



(16b) 



and by virtue of (TTTj) the boundary values of the species 
concentrations are 



_ T i(Xst) n 

°i,st — XstLf , 

a 

Ci, st = T ^(l-X*)G, 



(16c) 



It should be noted beforehand that the given model de- 
scribing, generally speaking, surface kinetics contains at 
least three variables, the solid state composition \ and 
two boundary values of the species concentrations Cf, 
Cf. So the system instability can be described using 
the classical notions of relaxation oscillations on a two- 
dimensional phase plane only at a rough approximation, 
as it has been already shown in our previous paper [Tll | . 



III. THE INSTABILITY DOMAIN 

A. The eigenvalue problem 

Now let us analyze in a rigorous way the linear stabil- 
ity of the system around the steady state described by 
expressions (|16[) . For this purpose the dynamics of small 
perturbations 

6Ci(t, z) oc exp {jt - piz} , 8x(t) oc exp {-yt} (17) 

in the species distribution and the composition of the 
crystal surface is considered. Here 7 is the instability in- 
crement and the parameters {pi} such that Rep; > 
characterize localization of the perturbations SCi (t, z) 
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near the crystal surface. Then the governing equa- 
tions (j9"j)- (fTTj) , and (fl~5|) are linearized with respect to 
perturbations ifTTj) in the vicinity of the stationary solu- 
tion p6p . The system of algebraic equations obtained in 
this way gives us the eigenvalue equation for the instabil- 



ity increment 7. This procedure is practically identical 
to that from Ref. [TlJ] . So here we skip the correspond- 
ing mathematical manipulations and write directly the 
desired eigenvalue equation in the final form 



r 2 
9 



-i + x(i-x) 



(CA)e 



(C/A)e^ 



(CA)e^ + 1 



(C/A)e"/' + 1 



(18) 



where following the notations of paper we have introduced the variable C > 0, the angle ip £ (— tt/2, 7t/2), and 
the parameter A > given by the expression 



with 



(19) 



such that 



and the diffusion flux takes the value 



7 = 



and 



Pi 



a 1 
D 2 T2{x) A 



i'0 



(20) 



(21) 



The quantity g stands for the dimensionless diffusion flux 
of species through the aqueous solution bulk towards the 
crystal surface 



g = ^D 1 D 2 T 1 (x)T2(x)G 

= VDiD 2T 2 exp {0(1 - 2 X ) + 0} G . (22) 

To find the boundary of the instability region in the 
space of system parameters we note that the eigenvalue 
equation (jTSJ) can be directly reduced to a fourth-order 
polynomial equation by multiplying it by both the de- 
nominators entering its right-hand side. The coefficient 
of the highest power term of this polynomial is a constant 
value. So the roots of equation (Tig)) cannot go to infinity 
and, thus, vary continuously as the system parameters 
change. The instability boundary separates the regions 
where the value of Re 7 has different signs and, therefore, 
meets the equality 

Re7 = 

converting, due to (|20| . into the condition ip — ±7r/4. 
Taking the latter into account and splitting equation (fl8|) 
into the real and imaginary parts we immediately get the 
conclusion that at the instability boundary the parameter 
C obeys the following equation 



(0 + $*i(CA) 



26 c 



(23) 



V29 c Cc 



)A* 2 (c c A)+(0-0)i* 2 (!) 



(24) 



where £ c is the solution of equation (|23|) and the functions 



T , , V2x(V2a; + l) 



1 



{V2x 



1) 2 + 1 

(25) 

as well as the critical value of the nonideality parameter 
depending on the crystal composition x 



e c (x) 



1 



2x(l-x) 



(26) 



have been introduced. In other words, at the instability 
boundary the general eigenvalue equation (|18p is reduced 
to (|23|) and if its solution £ c exists then formula (|24|) 
specifies the critical value of the species diffusion flux g c . 
Only one additional condition should be imposed; it is the 
requirement that the obtained value of g c be positive. 

Below we will confine our consideration to the case 
X = 0.5 only for which 9 C := 9 c o = 2. It due to, first, 
exactly this value of the solid composition x determines 
actually the boundaries, external and internal ones, of 
the instability regions to be analyzed. Second, as follows 
directly from expressions (|2"3"]) and by transforma- 

tions 



7old 



e, 



7,0 



e,. 



Did • 



7<:0 



(27) 



the case of x 0-5 is reduced immediately to the given 
one. Naturally the dependence of the system charac- 
teristics on the solid composition x endows the growth 
instabilities with nontrivial properties. In particular, in 
some sense "optimal" conditions of the instability onset 
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FIG. 2: Left-hand side (LHS) of equation (I23p as a function 
of the variable (. Curve 1 depicts this dependence when the 
system asymmetry cannot affect the instability onset crucially 
(4> < 6), curve 2 exhibits the case where the asymmetry effect 
is pronounced (<f) > 8). Curve 3 demonstrates the fact that 
the system asymmetry depresses the instability onset when 
d>> 9 and A < 1. 



can match the solid composition deviating substantially 
from x — 0-5, which in turn is able to cause a system 
instability with respect to spatially nonuniform pertur- 
bations. This question, however, is beyond the scope of 
the present paper. Third, in the mathematical expres- 
sions to be obtained below the quantity 9 c q will be kept 
on instead of being replaced by its numerical value, so 
using transformations (|2T[) the general expressions can 
be reconstructed immediately. 

The solution of the system lj2"3")) and implicitly de- 
termines the critical value g c (9, <f>, A) of the species dif- 
fusion flux. Thereby it describes the boundary of the 
instability region in the complete space of the system pa- 
rameters {g, 9, 4>, A}. Projecting this region onto various 
planes makes it possible to regard the instability bound- 
ary as some curve (or surface) dividing a given plane into 
two domains, where the instability can arise in principle 
for a given values of the corresponding parameters or can- 
not do it at all. Below in this Section we will consider 
in detail this projection onto the plane {9, </>} for a fixed 
value of the parameter A with the main attention paid 
to the limit A > 1. 



B. Two mechanisms of the instability 

Possible roots {Q} of equation (|2"3")) specify the eigen- 
values determining the critical value of the diffusion flux 
g c via expression (|24[) . The instability boundary is the 
locus where the potentials <j> and 9 take such values that 
the left-hand side of equation (J22J) gets its maximum at 
these roots. The solid nonidcality and the system asym- 
metry are responsible for the terms in this equation ex- 
hibiting different behavior. The term proportional to 9 
is the increasing function of £, whereas one proportional 



to 4> comprises increasing and decreasing branches. It is 
the mathematical reflection of different instability mech- 
anisms caused by the solid nonideality and the system 
asymmetry. 

The asymmetry effect becomes crucial when the left- 
hand side of equation (|2"3")l changes its behavior as a func- 
tion of £. It converts from a function monotonically in- 
creasing from to 29 when £ runs from to oo (curve 1 
in Fig. ^ to one possessing a maximum ^ m attained at a 
certain internal point < ( m < oo (curve 2 in Fig. [J). For 
< ( < ( m it grows from to \& m > 29 and then drops 
down to 29 on the interval Cm < C < 00 • The asymp- 
totics of the left-hand side of equation (f2"3")l as C — * 00 
demonstrates us directly that it is the case when 



A > 1 . 



(28) 



In fact, if inequality (|28p holds the asymptotics of the 
left-hand side of equation ([23]) is a decreasing function of 
£ and, thus, the point < £ m < oo does exist. Exactly 
in this case the instability can arise even the nonidcality 
potential is less then its threshold, 9 < 9 c q provided the 
maximum 9 m > 29 c q due to the effect of the system 
asymmetry. For the latter to be the case the inequality 
A > 1 is necessary as follows from condition (|2"8")l . For 
A < 1 the system asymmetry depresses the instability 
onset as it is illustrated in Fig. [3 by curve 3. 

Before passing to a detailed analysis of the instability 
domain we present a fairly simple way to construct the 
instability boundary of the plane {9, (f>} for a fixed value 
of A. It applies to the fact that the given system admits 
two scenarios of the instability onset. One caused by the 
solid nonideality matches the eigenvalues £e*^ — > oo with 
the diffusion flux g — > oo. In this case the solution of the 
general eigenvalue equation (|18[) can be written as 



for 



oo , 



(29) 



so the instability arises when the nonideality parameter 
exceeds its critical value, 9 > 9 c q, because 7 oc Q 2 e l2 ^ . 
The other is characterized by the bounded variations of 
the eigenvalues (, e as the diffusion flux goes to infinity. 
Under this condition we can analyze directly the eigen- 
value problem in the limit g — > 00 setting the left-hand 
side of equation (fT5|) equal to zero and, thus, reducing it 
actually to a quadratic equation. Omitting simple arith- 
metical manipulations the result is 



where 



1 



4(0 cO - 9) 



K ± ij K 2 



l69 c0 (9 cQ -9) , (30) 



K := A 



A- 



2 A 



(31) 
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FIG. 3: The instability boundary for several values of the 
parameter A including the limit value A = oo. The used 
criterion of instability is Re 7 > for g — > 00. 



For 9 > 9 c o one of these roots corresponds to unstable 
perturbations, nevertheless, the perturbations matching 
the eigenvalues given by expression (|29p are dominant 
due to large values of their increments. However, when 
the solid nonideality is not to high, i.e. 9 < 9 c q, the latter 
perturbations turn out to be stable and the growth in- 
stability is caused by the system asymmetry. Indeed, the 
instability boundary with ip — ±7r/4 meets the condition 



K > 



and 



89 c0 (9 c0 — 9) , 



(32) 



By virtue of (|32p such instability can arise when the pa- 
rameters 4> and A > 1 reflecting the system asymmetry 
meet the inequality 



> 



A 2 



A 2 - 1 



— (26> c0 — 9) + ^ \/ 9co(9co — 9) . 

When A < 1 the system asymmetry suppresses the insta- 
bility as noted above applying to Fig. [21 The curve 
on the plane {9, </>} specified by the dependence <j>t{6) is 
presented in Fig. [3] for several values of the parameter A. 
Roughly speaking B^ is the boundary of the instability 
domain for 9 < 9 c q. 

It should be underlined that the present analysis was 
based on the assumption that the instability has to arise 
for large values of the species diffusion flux if it can de- 
velop in principle for given values of the other system 
parameters. It is true when the growth instability is 
caused by the solid nonideality. However, for the insta- 
bility induced by the system asymmetry the situation is 
more intricate. Rigorously speaking, in the latter case 
at the real instability boundary B^, the species diffusion 
flux g takes a certain finite value g c < 00 and in a narrow 
boundary layer inside the instability region the diffusion 
flux must belong to a finite interval, g c < g < g£ < 00. 
Nevertheless, as will be seen below, this feature is valu- 
able only for A > 1. So as stems from expression (|3"3"]l 



FIG. 4: The structure of the instability region as a whole on 
the plane {9, c/>} for a fixed value of the parameter A 3> 1. 
It comprises five regions distinguishable in properties: two 
volumetric domains T>g and 2?o-i, one intermediate layer Ce 
between them and one boundary layer C4, whose thickness 
Wc > 1/A <§; 1, and, finally, a double criticality neighborhood 
C of the point {8 c0 , </> c o ~ c0 }. 



for A ^ 1 the boundary of the growth instability caused 
by the system asymmetry is approximated by the line 



4> = 26> c0 - 
at the leading order in 1/A. 



(34) 



IV. STRUCTURE OF THE INSTABILITY 
DOMAIN 



Based on the analysis of the eigenvalue equation (|23|) 
for A ^ 1 we can single out five characteristic regions of 
the system instability on the plane {9, 0} shown in Fig. [4] 
Let us consider them individually assuming A > 1 to 
hold. 



A. Instability domain Vg 

The volumetric domain T>g matches actually the 
growth instability studied in our previous paper [11| . It 
is bounded by the vertical line Bg = {9, <f) : 9 = 9 c0 = 2}, 
by the layer Cg, and the f9-axis. The layer Cg is a certain 
neighborhood of the line (j> = 29 c0 — 9 whose thickness 
is about W c > 1/A. In domain Vg condition (12"5)) is 
strongly violated, i.e. 



A 2 > 1 



(35) 



So the left-hand side of equation ([23)1 is monotonically 
increasing function of C- Besides, for any point of the do- 
main T>g the distance between it and the line </> = 29 c q — 9, 
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i.e. actually between it and the layer Cg, can be re- 
garded as a large value in comparison with the quantity 
1/A. The latter statement, as can be shown directly, 
causes the solution of equation (f2"3"| to meet the inequal- 
ity Cc 3* 1. Thereby the former term on the left-hand 
side of equation (123|) can be taken in the limit (A — > oo. 
Also the corresponding term in expression (|24|) is ignor- 
able. Under these conditions equation ([2"3")) is reduced 
to a quadratic equation with respect to £, yielding us 
immediately its solution in the form 



Cc 



V2 



F{rg) 



(36) 



Here, by definition, the function F(rg) is determined by 
the expression 



F(x) 



1 



2(1 - x) 



2x - 1 + v/1 + 4z(l - x) 



its argument is 



26*co — $ 



1 



(37) 



(38) 



The inequality 9 c q < 9 is assumed to hold, thus, < 
rg < 1. Then the critical value <jf c of the dimensionless 
diffusion flux is 



9c{T> e } 



A 2 c0 F\rg)[F{r e ) + l] 



re 



(39) 



by virtue of (124)) . 

Near the threshold of the nonideality coefficient, i.e. in 
the vicinity of the boundary Bg, where 



< 6- 



«0- 



(40) 



one has 1 — rg <C 1. In this case function (|37j) is ap- 
proximated as F(rg) ps 1/(1 — rg) and expression (|3^|) is 
reduced to 



5c{x> e |B 9 } 



?c0 1 



(41) 



Whence it follows, in particular, that the critical value 
of diffusion flux diverges as (9 — 9 c o)~ 3 for 9 — > 9 c q + 0, 
being in agreement with the results of paper 

[HI- 

When the analyzed point {9, </>} is located in a close 
proximity to the layer Cg, i.e. for 



< 29 c 



(42) 



one has rj < 1 and F(rg) « 2rg. In this case for- 
mula (|39p is simplified as 



9c{T> e \Ce} 



4A 2 cO (29 cQ - 9 



(43) 



We remind that in expression (|43p the difference (20 c — 
9 — cj>) cannot become too small according to inequal- 
ity (|42|) . The behavior of the critical diffusion flux for 
points coming close to the line <j) = 29 c0 — 9 is considered 
below. 



B. Instability domain X>o-i 

The other volumetric domain 2?o-i of system instability 
is formally the half-plane bounded from below by the 
layer composition Cg [J Cs, i-e. by a neighborhood of the 
line 4> = 29 c0 - 9 with thickness W c > 1/A (Fig. gj. It 
comprises all the points {9, 77} such that 



<t> - 29 cQ > — 



(44) 



In this region the solution C c of equation (|23|) turns out 
to be much less than unity, ( c C 1. As can be verified 
directly the latter inequality enables us to ignore both 
the second terms on the left-hand side of equation (T2"3"]) 
and inside the square brackets in equation ([24]) . The ap- 
pearance of these terms is due to the channel 0-2 of the 
precipitation reactions (fT|) . Therefore in the domain 2?o-i 
the contribution of the channel 0-2 is of minor impor- 
tance and the growth instability is caused by the channel 
0-1 individually. 

Using this simplification the eigenvalue equation (|2"3"1) 
again can be reduced to a quadratic equation with the 
solution 

Cc = ^*W), (45) 

where, by definition, the argument roi is the value 

2# c o 



(46) 



and meets the inequality 1 — r i ^> 1/A. 

The corresponding expression for the critical value of 
the species diffusion flux is 

1 



{v - 1 } = ^F 2 (r 01 )[F(roi) + l] 



(47) 



In particular, near the domain boundary Cg [J C^, i.e. for 
1 



< 9 + - 29 c0 « 1 , 



(48) 



where 1 — roi <C 1 and the function F(roi) ~ 1/(1 — roi) 
expression (1471) converts into 



9c{T>o-i\C e \JC<i,} 



'c0 



A 2 (6>- 



29. 



cO) 



(49) 



It should be pointed out that expression Ij49|) does not 
describe a real singularity in the diffusion flux threshold. 
In fact, the difference 



29 co 



(50) 



is bounded from below in the domain 2?o-i 1 namely, a 3> 
1/A and the limit a — > +0 cannot be implemented in it. 

By virtue of expression (|47| in the domain T>o-i the 
diffusion flux threshold g c practically does not depend 
on the particular value of the difference (9 — 4>) because 
of the minor effect of channel 0-2. The situation changes 
dramatically when the analyzed point {9, <f\ enters the 
boundary of this domain, being the subject of the follow- 
ing subsections. 
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FIG. 5: The intermediate layer Ce separating the instability 
domains T>g and X>o-i. Zoomed-in view. 



C. Intermediate layer Co 

The instability domains T>g and 2?o-i are joined to each 
other via the layer Ce whose points are located near the 
line 9 + <j) = 29 c q and meet the inequality 9 > 4>. So 
the left-hand side of the eigenvalue equation ([25]) is a 
monotonously increasing function of £ and the solution 
of this equation £ c decreases as the potential 9 increases 
for a fixed value of <p. The results obtained in the two 
previous subsections show us that the quantity £ c , first, 
drops from very large values up to £ c > 1 as the analyzed 
point {9, <f>} goes from the instability boundary Be to the 
layer Ce- Then, just after the point crossing the layer Ce 
the quantity Ce takes values about C, c < 1 and drops down 
to zero as the analyzed point goes away from it. In fact, 
on one hand, by virtue of (|36[) when the analyzed point 
{9, cf>} tend to the layer Ce on the side of the domain T>e 
and inequality (jl2")l holds we have 



Cc 



\/2A(20 c o -9- 



On the other hand, for the point {9, (/>} located near the 
layer Ce on the side of the domain 2?o-i where the in- 
equality (|48|) hold the solution £ c of the eigenvalue equa- 
tion (|2"3")l is approximated as 



A> 1. Namely, in this subsection we consider the region 
for which 



29 c0 \ « 1 



thereby the two inequalities 



C C A > 1 and 



< 1. 



(51) 



(52) 



hold simultaneously. This region comprises the layer Ce 
as well as the neighboring parts of the domains T>e and 
2?o-i- So the expression for the diffusion flux threshold 
valid in it really specifies the crossover between the do- 
mains T>e and X>o-i • 

Under condition (p)2")) the former term on the left-hand 
side of equation (|23p can be approximated by the asymp- 
totics of the function ^i(x) for x — > oo, whereas the lat- 
ter one matches the limit x — > 0. Therefore in this case 
equation (|23p can be rewritten as 

(0 + <£)-- {9-<f>)C= V2Aa, (53) 



where a is given by expression (|5tJ|) . The solution of ([5 
is of the form 



Cc 



1 



V2(9- 



- 4?) - Act 



Then the substitution of J54]) into ((25) yields 



5c{£(,} 



20 cO AC c 2 



VA^2 



2 ) 



(54) 



(55) 



where the function ^2(2;) has been also approximated 
using the appropriate asymptotics. As it must be ex- 
pression (f5"B"|) converts into expressions (|43|) and (TJ5]) for 
=FAct 3> 1, respectively. 

Figure [6] illustrates the obtained crossover of the dif- 
fusion flux threshold. In this Figure the critical value g c 
of the dimcnsionless diffusion flux is shown vs actually 
the nonideality parameter 9 for fixed parameters <fi = 1 
and A = 50. In other words, it visualizes g c (9) for the 
analyzed point {9, <f\ moving along the line te shown in 

Fig.m 



D. Boundary layer C 4 



V29 c0 
U A(9 + <j>- 29 c0 ) 

by virtue of ([45]) . The "boundaries" of the layer Ce meet 
the estimate A\9 + <fi — 29 c q\ > 1, which justifies the 
statement mentioned above. So inside the layer Ce the 
quantity £ c has to change in the interval 1 < Cc ^5 1- 

The expression obtained below for the critical value of 
the diffusion flux g c is valid, however, for a wider region 
than the layer Ce itself due to the adopted assumption 



When the potential 9 is less than the threshold 9 c q 
the solid nonideality cannot individually cause the sys- 
tem instability. In this case only the cumulative effect 
of the system nonideality and asymmetry gives rise to 
the growth instability or even the system asymmetry it- 
self does when the potential is high enough. So for 
9 < 9 c q the instability domain 2?o-i borders with the re- 
gion of the stable crystal growth via the boundary layer 
C<jj (Fig. [7]). Let us consider its properties in detail. As 
for the layer Ce analyzed in the previous subsection the 



9 



100 




FIG. 6: The critical value of the dimensionless diffusion flux g c 
vs the nonideality parameter 9 near the intermediate layer Ce . 
The plot is based on the general equations (|23|l and (|24|l using 
the parameters shown in inset. The straight lines visualize 
the formal asymptotics (|43[) and (|49[) whereas the dotted line 
corresponds to (f55l) . 




FIG. 7: The boundary layer Ce separating the instability 
domain 2?o-i and the region of the stable crystal growth. 
Zoomed-in view with an additional magnifying lens showing 
the finite structure of the instability boundary B$. 



given layer matches the root Q of equation (j2li)) of order 
unity, £ c ~ 1. However in this case by virtue of condi- 
tion (|28p the potential cj> should exceed the nonideality 
parameter, <f> > 9, for the growth instability to arise. As a 
result the two terms entering the left-hand side of equa- 
tion (|23| have opposite signs and the functions $>i(x), 
^2(2:) should be approximated to the next order in the 
corresponding small parameters in comparison with the 
case of the layer Co . Namely, using inequalities (|52p ex- 
pressions (|23|) and (|24|) are reduced to the equation 



(4 + 0) 



6)C = V2Aa 



(56) 
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FIG. 8: Illustration of the mechanism responsible for the com- 
plex behavior of the critical diffusion flux in the boundary 
layer Cp (upper fragment) and the resulting i? c (#)-dependence 
(lower fragments). The right-hand fragment depicts this 
dependence in zoom, making it evident that in the region 
(#ci, 0c2) the growth rate should belong to a bounded interval 
for the instability to arise. In plotting the potential differ- 
ence a as a formal function of the variable ( determined by 
equation (|56[) and the function $(C) (see equation l|57p) vs 
the variable ( the values 9 = 1.5, (j> = 2.5, and A = 10 were 
used as well as the original functions 9i(x) and 3*2 (z) (ex- 
pression (|25p ) rather then their approximations were applied 
to take into account not only the leading terms but also all 
the other small contributions. 



and the expression for the diffusion flux threshold 



*(0 



(57) 



where the function 3>(C) is introduced by the formula 



*(0 



V2 
A 



+ > 



C 2 



9)C 



o 



1 

A~2 



(58) 



It should be pointed out that expression (|56|) does not 
contain a term of order A -1 and the other term of order 
A~ 2 is not written explicitly because its effect is reduced 
only to a small constant contribution to the value of a. 



10 



To explain the resulting dependence of the diffusion 
flux threshold g c on the potential difference a let us refer 
to Fig. [H It illustrates the value of a treated as a formal 
function of £ that is determined by equation (|56[) and has 
a minimum a m meeting the estimate 



(59) 



and attained at 




(60) 



It matches the critical value of the species diffusion flux 

e c0 A 2 (<t> + e\ 

9m = —x— (61) 



written in the leading order of 1/A. 

If the potential difference a = + 9 — 29 c o is less than 
a m there is no solution of equation ([56]) (i.e. of equa- 
tion (|23p ) and the crystal growth is stable. When a > a m 
equation (|56p admits two solutions written again in the 
leading order of 1/ A as 



C = V2(0 + 9) A<j+ v/A 2 a 2 - 2(0 2 - 9 2 ) 



and 



V2(0 - 9) 



Aa+ ^A 2 ct 2 - 2(0 2 - 6» 2 ) 



(62) 



(63) 



Solution (|62|) matches the decreasing branch of the de- 
pendence cr(C) (Fig. [8]) and describes the lower boundary 
of the diffusion flux threshold g c {&) obeying the estimate 



1 



28^ 



g c (a) \ 2U ■ "7 i cr„ 

' (64) 

This branch actually specifies the minimal value 
9 c i(<p, A) of the nonideality parameter 9 necessary for 
the growth instability to arise for given values of the pa- 
rameters <fi and A, namely, by virtue of (I5ii 



29 r 



2V2 



(65) 



As it must be expression (|64|) converts into expres- 
sion P9"]) for Act ^> 1 describing the behavior of the 
diffusion flux threshold in the instability domain X>o-i 
near the boundary layer 

Solution (|63p describes the upper boundary of the in- 
stability region g+ (a) which, however, exists only within 
a rather narrow interval of the potential difference er, i.e. 
when a m < a < cr m (Fig. [5]). The parameter cr m and the 
corresponding value Cm niatch the point where the func- 
tion ^(C) changes its sign passing through zero. Exactly 



at this point the upper branch 5+ (a) of the diffusion flux 
threshold goes to infinity and for ct > cr+, i.e. for C > C m 
it does not exist. In this case the values of the species 
diffusion flux corresponding to the instability onset are 
bounded only from below by the threshold g c {&)- Ac- 
cording to expression ([5"6")) - (j5"7| the difference between 
Cm and Cm is a value of the first order in the parameter 
1/A, namely, 



Cm Cn 



y/29 
A(</>-0) 



(66) 



and as a result the corresponding difference of the pa- 
rameters is 



5a n 



V2 9 2 



AV<^ 2 -f 2 



(67) 



The obtained expression demonstrates the fact that this 
difference specifying the thickness of the region where 
the diffusion flux threshold exhibits complex behavior is 
extremely narrow (Fig. [H]). It is of the third order in 
the small parameter 1/A and can be ignored at all. In 
this case only the first term in expansion (f58|) should be 
taken into account, thus, only branch (|59p exists. So by 
virtue of ([r3| the diffusion flux threshold in the layer L$ 
as well as in its small neighborhood meeting the interval 
<7 m < cr <C 1 is approximated by the expression 



9c(<j) 



160 c 3 o 



(68) 



showing some formal singularity when a — > a m + 0. 

Finalizing this subsection let us discuss the behavior 
of the instability boundary B$ depending on the param- 
eter A including its relatively small values. It should be 
reminded that previously we considered two curves on 
the plane {9, <fi}, the instability boundary = {(j> c (9)} 
itself and the curve B~f = {<p^(9)}. The former singles 
out the points on this plane where the crystal growth can 
become unstable for some values of the species diffusion 
flux. The latter is the boundary of the growth instability 
for vary large values of the diffusion flux, g — > 00. 

It can be demonstrated analyzing directly the general 
eigenvalue equation (f2"3"]l and the expression (f!M)) for the 
critical diffusion flux g c that the terminal points of the 
curves B^ and B^ at 9 = and 9 — 9 c0 coincide with 
each other for a given value of A. So by virtue of expres- 
sion (|3"3"]) their coordinates are specified by the following 
expressions 



= 29, 



A 2 + V2A + 1 
A 2 + l 



OcO 



7c0" 



1 



for 9 = 
for 8 = 8 c0 . 



(69) 
(70) 



As it must the ^-coordinates of both the points have a 
singularity as A — > 1 because in this limit the growth is 
stable for 9 < 9 c q. 
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FIG. 9: The difference in the instability boundaries B^ and 
B4, for several values of the parameter A. 




FIG. 10: Structure of the instability region in a close proxim- 
ity to the double critical point {6 c o, 4>co}- 



For the intermediate points < 9 < 9 c q the curves B<j,, 
B~f deviate from each other. To evaluate this difference 
Fig. [9] plots the difference 0+ — <j) c vs the potential 9 for 
several value of A. As seen in Fig.[9]the curves B^ and B~£ 
practically coincide with each other except for the values 
of A coming too close to its threshold A = 1. Thereby 
expression (|33l) gives a fairly fine approximation of the 
instability boundary B$ for such values of A. 



E. Double critical point and its neighborhood C 

The boundaries Bg and B^ of the instability region 
meet at the point {# c o,</> c o} that can be referred to as 
a double critical point because its coordinates are the 
threshold of the nonideality parameter and the thresh- 
old of the asymmetry potential exceeding which the sys- 
tem asymmetry changes the instability property substan- 
tially. The latter implies the fact that the asymmetry 
causes the instability onset in the system being stable 
before the potential exceeds the threshold, <f> > (j> C Q, and 



4> C Q is the minimal value possessing this property among 
all the possible values of the solid composition \ and the 
nonideality parameter 9. Therefore in calculating the 
value of 4> c o we can set 9 = 9 c q . 

The critical region C is a certain neighborhood of the 
point {6> c o, 4>co\ where the layers Le and L<f, overlap with 
each other. So it should exhibit some crossover between 
the properties of these layers. According to the results to 
be obtained in the region C the potential difference <f>— 9 is 
rather small so not only the inequality £ C A 3> 1 but also 
Cc/A ^> 1. Keeping in mind the general condition l|28p 
necessary for the system asymmetry to affect essentially 
the instability onset we describe the region C with two 
small parameters u«l and v <C 1 introduces as follows 



4>. 



A 2 



-(! + «), 



?c0 



(71) 



2V2A 2 



Then for the variable £ :— A/£ <C 1 regarded as a small 
value the eigenvalue equation (|2"3")l is reduced to 



-xu + X 



(72) 



and expression 
the form 



for the diffusion flux threshold takes 



9c{c\ 



l (xV2- 



(73) 



As it must, when u < the instability boundary is spec- 
ified by the equality v — (9 = 9 c q) and the diffusion 
flux threshold g c — > oo as v — > +0. For u > the system 
changes the behavior. 

The eigenvalue equation (|72p relating the variables u 
and v at the point x = y/u/3 where its right-hand side 
attains the minimum specifies the instability boundary 
Bs, namely, 



3%/^ 



.3/2 



(74a) 



or returning to the variables 9 and 



7<:0 



6a/3 



3/2 



(74b) 



As should be expected, at the boundary B$ the diffusion 
flux threshold takes a finite value equal to 



Sc{C|£v} 



6a/6A 



3/2 



(75) 



Naturally, the diffusion flux threshold g c diverges as the 
asymmetry potential 4> — * ^c0 + 0. 

Near the boundary B^ the values of the diffusion flux 
causing the instability onset are bounded from below and 
above. The locus B~t where the upper boundary goes 
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FIG. 11: Perturbation of the species distribution in the aque- 
ous solution induced by variations in the surface concentration 
C'f on time scales about r. Schematic illustration. 



to infinity is specified by the singularity point of func 
tion (f73|) . i.e. x = u/v2. This value via equality lf?2] 
gives us the relationship between the potentials 9 and c, 
at the curve Bt 



1 2 

V = =u 

V2 



A-0 



'<■() 
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(76a) 



(76b) 



The expressions obtained here hold for u,v <C 1, so the 
characteristic size of the region of double criticality is 
about R c ~ 1/A 2 . 



V. REGIMES OF INSTABILITY DYNAMICS 

The present section is devoted to a qualitative analysis 
of the system dynamics. For the sake of simplicity we 
ignore difference in the species kinetic coefficients setting 
D\ = Z?2 = D and v\ = V2- 

At first, let us consider perturbations of the species 
distribution 5Ci(z,t) induced by small variations Sx{t) 
in the surface composition on time scales about r. Ac- 
tually 1/r is the perturbation increment analyzed in the 
previous Section. Change in the surface composition x(t) 
affects directly the species attachment rate caused by the 
growth process, which, in turn, gives rise to variations in 
the species concentration near the crystal surface SC-. 
These boundary variations in the species concentration 
spread into the solution bulk, which is responsible for the 
formation of spatial perturbations in the species distribu- 
tion schematically shown in Fig. [TTJ The characteristic 
spatial scale of these perturbations can be estimated as 
h T ~ {Dt) 1 ' 2 . 

Within a qualitative approximation mass conservation 
for such perturbations reads 



h T 5C! 



S[n(C!, x )] 



(77) 



or, by virtue of (|T2|) . 



!^5C* ~ -5Cf 

T Ti 



aCi 



where the quantities (for i = 1, 2) 

dlnrj(x) 



Wi(Xst) 



dx 



have been introduced and by virtue of {T3 



(78) 



(79) 



(80) 



Expression (|75|) enables us to single out two limit cases. 
The first one which will be referred to as the growth 
regime of constant growth rate matches rather slow vari- 
ations of the crystal composition x and the species con- 
centration Ci, namely, the condition r 3> Ti(h T /a) or, 
what is the same, 



r > 



Dt? 



In this case (|75|) yields 

sc? 

and, thus, via ([77]) 



-Cl st uJiSx 



S[n(C!, x )} 



Dt 



1/2 



< r 



2,st 



(81) 



(82) 



(83) 



Thereby for slow variations of the crystal composition 
X the induced perturbations in the species distribution 
Ci{z,t) are in quasiequilibrium. In other words, the 
boundary value C? of the species concentration changes 
in time with x m such a manner that the boundary 
value of the diffusion flux, the species attachment rate 
ri , be practically equal to the inflow of the corresponding 
species at distant points. In particular, exactly such vari- 
ations are described by expression ()82|1 being lineariza- 
tion of the condition 



aC[ 

n{x) 



const . 



(84) 



The second limit case, which will be called the growth 
regime of constant surface concentration is related to 
rather fast variations in the crystal composition when 
their time scale r meets the inequality r <^ Ti(h T /a) or 



r < 



Dt? 



(85) 



In this case the induced variations in the surface con- 
centration Cf of species i are rather small in comparison 
with that could be expected in the previous limit case, 



6C? 



Dt? 



1/2 



■0HCl at 8x<u i Cf 8t S X . (86) 
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Therefore, at the first approximation the fast dynamics 
of the species distribution and the crystal composition 
meets the equalities 



Cf 



const and <5r, 



iCfd 



1 



n(x) 



(87) 



In this consideration the variations of the crystal com- 
position x(t) were treated to be given beforehand. In 
order to draw some conclusions about the growth dy- 
namics as a behavior of an autonomous system it is nec- 
essary to discuss how the induced variations of the at- 
tachment rates r\ and r-i affect, in their turn, the crystal 
composition \- This effect is described by the governing 
equation ([To"} . 

As was discussed in the previous Section, when the 
nonideality potential exceeds the critical value, 9 > 9 c q 
the perturbation increment 1/r — > oo as the species dif- 
fusion flux goes to infinity also. So it is natural to ex- 
pect that for the developed instability the regime of con- 
stant surface concentration takes place with respect to 
both the species components. Then keeping in mind ex- 
pressions ([57} and applying to equation (TT5j) governing 
the dynamics of crystal composition we can draw the 
velocity field of the system motion on the phase plane 
{Cf/Cf,%} as shown in Fig. [121 The curve 



Cf 
C 9 S 



X T i{x) 



X 



(i-x)t 2 ( x ) (1-x) 



divides this phase plane into parts with the opposite di- 
rections of the velocity field. In obtaining ([88} expres- 
sions ([13} have been used. As it should be the stationary 
values of the species concentrations Cf ai and the crystal 
composition % s t (see expressions dM}) meet equality (1551 . 
Figure [12] clearly demonstrates us that under such condi- 
tions its increasing branches are stable whereas a decreas- 
ing branch (if it exists) is unstable. So the limit circle 
at a rough approximation should have the form shown 
in Fig. 1121 Fjxactly this limit was analyzed in our previ- 
ous paper [Tlj and corresponds to the domain Dg of the 
instability region. 

If the nonideality parameter 9 is less then the criti- 
cal value, 9 < 9 c q — 2 the solid nonideality cannot itself 
induce the growth instability. In this case the instabil- 
ity development is governed by the system asymmetry, 
which is reflects in properties of the instability domain 
2?o-i- I n particular, for the system with such parameters 
only the channel 0-1 of the precipitation reactions (fT} 
plays an active role, the channel 0-2 is characterized by 
the equilibrium value of the species diffusion flux at the 
crystal surface. In this case it is quite natural to assume 
that the perturbation increment 1/r meets the inequality 



Dt 2 Dtx 

— S- < T < 



(89) 



Therefore, on one hand, with respect to species 2 such 
a process can be classified within the regime of constant 
growth rate. On the other hand, with respect to species 1 



is 



I. solid nonideality effect 

II. system asymmetry effect 
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FIG. 12: Phase planes demonstrating the mechanism of the 
instability onset caused by the solid nonideality (I) and the 
system asymmetry (II). 



the regime of constant surface concentration takes place. 
Actually it is the case for the points of the domain 2?o-i • 
To describe the corresponding dynamics of the crystal 
composition \ we can fix the surface concentration Cf 
and set the species attachment rate r% = (1 — Xst)G. 
Then we draw a similar velocity field of the system mo- 
tion on the phase space {Cf , x} shown again in the same 
Fig. [12] Its pattern is identical to one discussed above 
except for the fact that the y-axis of this phase plane 
has now another meaning, it presents the surface con- 
centration of species 1. As follows from equation (fn 
and expressions ([T3} the curve 



(i-x) 



n{x)XstG oc 



(i-x) 



-{e+4>)x 



(90) 



separates the regions on the phase plane {Cf,x} with 
the opposite directions of the velocity field. This curve 
looks like the previous one ([88} within the replacement 
29 — ► 9 + (f>. So again the instability condition for 
the potentials of the species interactions take the form 
9 + <p > 29 C Q, being in agreement with the results ob- 
tained before. As previously the increasing branches of 
curve (|90p are stable whereas the decreasing one is unsta- 
ble and the system transition between them as well as the 
transition from the unstable stationary point {Cf gt , Xst} 
to one of them proceeds within the regime of constant 
surface concentration with respect to species 1. The 
rough approximation of the limit circle again has the 
same form. 

In the part of the domain I\>-i where 9 > 9 c o both of 
the instability scenarios can be implemented. So depend- 
ing on the species diffusion flux either the phase plane 
{Cf , x} or the plane {Cf /Cf , x} can gi ye an appropriate 
representation of the system dynamics. 
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VI. NONLINEAR DYNAMICS OF SYSTEM 
INSTABILITY. DOMAIN V$ 

This section presents numerical results for the system 
dynamics when the growth instability arise in a subdo- 
main of the domain £>o-i, where the nonideality pa- 
rameter 8 is less then its threshold, i.e. 6 < 9 c q. So it is 
the the system asymmetry that causes the instability. 

To model numerically the system dynamics the gov- 
erning equations (|9j)- (fTTj) , and (p~5|) were converted into 
dimensionless form. Namely, first, the time t and the 
spatial coordinate z are measured in units 



(91) 



respectively, i.e. the dimensionless time and spatial co- 
ordinates are introduced as t ncw — t Q id/T* and z ncw — 
z \d/z* . Second, the species concentrations and the diffu- 
sion flux are replaces with their dimensionless analogies, 
Cj, new = Ci )0 id/C* and G^ ncw = G^oid/G*, where 



C* 



1 



G* 



1 



(92) 



\/DiD 2 T g a 

In this way the original model is rewritten in the form 

d 2 a 



dt 

dx = 
dt ' 

with equation 
tion at z = 



8z 2 



Qi{x)C{-XQ2{x)C s 2 



(93) 
(94) 



being subject to the boundary condi 

da 



dz 



ft(x)C? 



(95) 



2=0 



and the condition at distant points, i.e. at the formal 
external boundary L nev! = L \&jz* \ 



Gi 



eg 

dz 



Here the dimensionless species diffusivities are 



m = — = 
k 2 




(96) 



(97) 



and the dimensionless rates of the atom attachment to 
the growing crystal are 



Qi{x) = ^iexp{<?(>x - 6(1 - x)} , 
Q2(x) = ^exp{-0(l - x) ~ Ox)} 



(98) 



with 



Kl = 



1 

^2 



1/2 



exp 



(99) 
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x = 0.5; G = 10 
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FIG. 13: The time dependent component SCi(t,z) of the 
species 1 distribution in the aqueous solution bulk near the 
crystal surface, z = 0, for several time slices within one pe- 
riod of the oscillations. Result of numerical simulation. The 
shown time origin t = is placed at an arbitrary chosen point 
that corresponds to the instability becoming well developed. 



It should be noted that the previously used parameter 
is related to the introduced kinetic coefficients as 



. 1/2 



(100) 



So the ratio (x\/ ' yc-^e^ is actually the main small param- 
eter of the given model because for aqueous solutions the 
relationship D\ ~ D 2 is typically fulfilled. 

The system of equations ([93 |) ~ (j96|) was solved numeri- 
cally using the Crank-Nicholson scheme for the diffusion 
equation and the midpoint method for equation (I94[) . 
To exemplify the basic characteristics of the instability 
dynamics in the region the system parameter were 
set equal to 6 — 1.5, <\> = 3.5, and = 10 as well as 
K\ = /ta = 1. Then expression (|100[) gave us the values 
of x\ and x 2 ■ The time and spatial steps in the simula- 
tion routine were 0.01, decreasing the steps twice did not 
affect the obtained results. The time variations in the 
species distribution induced by the developed instability 
turned out to be located near the crystal boundary within 
a layer of thickness about 15-20 spatial units. So the ex- 
ternal boundary of the system was placed at L = 100, 
where the species concentrations Cf° and C£° were fixed 
in such a way that the total diffusion flux and the solid 
composition take the values G Bt = Gf + Gf = 10 and 
X st = 0.5 under the steady state conditions. The total 
simulation time was 10000 time units. 

Below we will present the obtained results. Figure [T3l 
visualizes evolution of the species distribution in the 
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time f, in units of x* solid composition, x 



FIG. 14: The dynamics of solid composition x(t) an d the surface species concentrations Cf(t) (left column) and the corre- 
sponding phase portraits on the planes {X)C?} (right column). Result of numerical simulation. The shown time origin t = 
is placed at an arbitrary chosen point that corresponds to the instability becoming well developed and steady state. 



aqueous solution bulk near the crystal surface. Only the 
distribution of species 1 is shown because it, first, exem- 
plifies similar effects for species 2 also and, second, plays 
the leading role in the instability onset. To elucidate the 
dynamics of the species distribution the time dependent 
component 8C\(t, z) is singled out from the total distri- 
bution function 

d(M) = SC x (t,z) + (C{) + ^-z 

and depicted in Fig. 1131 The other terms in this ex- 
pression are the steady state components of the species 
distribution. As seen in this figure the time variations of 
species distribution are located near the crystal surface 
z = in its neighborhood of thickness about Lq ~ 15 for 
the chosen system parameters. So the size of the system 
L = 100 used in the numerical simulations is fairly large 
to enable one to regard the external boundary z = L as 
infinitely distant points. In any case in numerical simu- 
lations the size of the system should be specified that the 
inequality Lq < L to hold. 

Figure [13] demonstrates us the fact that a simple model 
of the boundary layer similar to the one shown in Fig. [11] 
can be used only for a qualitative analysis. The actual 
spatial form of 8Ci(t,z) can possess a remarkable ex- 
tremum attained at a certain internal point of the crystal 
neighborhood, which must be taken into account in con- 
structing an appropriate boundary layer approximation. 

Nevertheless, in spite of a rather rough model for the 



boundary layer used in Sec. |V] the instability scenarios 
described there is justified by the results of numerical 
simulation. The found dynamics of the solid composition 
X(t) and the surface species concentrations C-(t) exhibit 
relaxation oscillations with clearly visible fast and slow 
stages of system motion (Fig. [TJ]). So the results ob- 
tained for the given set of parameters do describe an es- 
sentially nonlinear regime of the growth instability. The 
phase portrait of the system oscillations on the plane 
{x, Cf } demonstrates us the fact that the regime of con- 
stant growth rate really takes place with respect to the 
species 2. Indeed the image of the oscillation limit cir- 
cle on this phase plane is located in the vicinity of the 
curve A^(x) obtained by setting the right-hand side of 
the boundary condition (|95p equal to the diffusion flux 
of species 2 under the stationary conditions, i.e. 

Q2(x)C^ (1 - Xst)G 
and thus 

N 2 (X)= (1 ~ Xst)Ge ^ .exp{-(0-g) x } . (101) 

With respect to species 1 the regime of constant sur- 
face concentration could be expected to be the case. The 
image of the oscillation limit circle on the plane {x, Cf } 
(Fig. IT4|) justifies this expectation at least within semi- 
quantitative consideration. Figure [HI depicts the ob- 
tained limit circle together with the nullcline N±(x) con- 
structed by setting the right-hand side of the governing 
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-^12 (x) determined by the expression 



0.2 0.4 0.6 0.8 
solid composition, % 



1.0 



FIG. 15: The phase portrait of the system oscillations on the 
plane {%, (Cf/Cf)}. Results of numerical simulation. 



equation (|94p equal to zero, fixing the surface concentra- 
tion Cf and assuming the attachment rate <?2(x)Cf of 
species 2 to meet the regime of constant diffusion flux. 
In this the expression 



(l- Xs t)Ge e X exp{-(0 + 0) X } 



(i-x) 



(102) 



has been constructed. As seen, here the fragments of the 
limit circle matching the fast motion deviate substan- 
tially from the decreasing branch of the nullcline Ni(x) 
and the fragments of slow motion go near its increas- 
ing branches. So, roughly speaking, it is the character- 
istics of the nullcline Ni(x) that specify the amplitudes 
of time variations in the solid composition and surface 
species concentrations for the developed growth insta- 
bility. However, the obtained limit circle also deviates 
remarkably from a simple form constructed in Fig. (reap- 
plying directly to the notions of the standard relaxation 
oscillations. The matter is that the system under con- 
sideration is really not reduced to a two-variable model 
implying actually the too simple boundary layer approx- 
imation shown in Fig. [TJJ to hold. So the dynamics of 
the surface concentration Cf (t) of species 1 contains the 
fragments of slow motion as well as that of fast motion 
(Fig. . The latter ones actually force the fast motion 
branches of the limit circle to deviate remarkable from 
horizontal lines on the plane {x,Cf}. This effect was 
also observed for the growth instability caused by the 



solid nonideality [11 1 



Finalizing the present Section we underline once more 
that there is a widely used approach to constructing the 
limit circle of oscillations in such system, i.e. the "bound- 
ary reaction - diffusion" systems treating the governing 
equation (|94|) (or its original version (p~5|) ) for the solid 
composition in a too simple way. It sets the right-hand 
side of this equation equal to zero and relates the sys- 
tem portrait on the plane {x, (Cf/Cf )} to the nullcline 



Ni2(x) 



(i-x) 



Q2(x) 
Qi(x) 



X exp{fl(l-2 X )} 

(i-x) 



(103) 



For the growth instability caused the solid nonideality 
the nullcline iVi2(x) possesses a decreasing branch being 
unstable (see, e.g. Fig. [T^J) - In this case the limit circle 
constructed following the classical ideas of the standard 
relaxation oscillations is justified at least within a quasi- 
qualitative analysis [Tlj . However, if the growth insta- 
bility is induced by the system asymmetry, such an ap- 
proach is not justified at all, the corresponding nullcline 
■^12 (x) is a monotonous curve and the system portrait 
on the plane {Xi(Ci/Cf)} is just located in its vicinity 
(Fig. US]). 



VII. CONCLUSION 

We have analyzed the oscillatory zoning, i.e. the self- 
organization phenomenon arising during crystallization 
of multi-component solid from aqueous solution. It man- 
ifests itself in self-formation of quasi-periodic spatial pat- 
tern of solid composition from the core of a crystallite to 
its rim. 

Keeping in mind systems like (Ba,Sr)S04 we have pro- 
posed a model for the growth of ternary-component solid 
from aqueous solution. The crystallization process com- 
prises passive diffusion of species towards the crystal sur- 
face through the aqueous solution bulk, their adsorption 
at the crystal surface, and incorporation into the crys- 
talline lattice at the surface atomic steps. The latter pro- 
cess is assumed to limit the crystal growth, so the species 
adsorption-desorption at the crystal surface is described 
within the quasi-equilibrium approximation. Due to a 
very low rate of crystallization from aqueous solutions 
the growth dynamics is simulated using the boundary- 
reaction-diffusion model for the species distribution in 
the aqueous solution bulk. 

The proposed model for the growth process takes into 
account the solid nonideality as well as the system asym- 
metry, with the latter being the characteristic feature of 
systems for which oscillatory zoning was reproduced in 
laboratory also. It has been demonstrated that the sys- 
tem asymmetry can cause the growth instability in the 
case when the solid nonideality is low, i.e. the nonideality 
parameter is less than its threshold, 9 < 9 c q, or even if 
the solid solution is ideal, 9 = 0. Using the linear stabil- 
ity analysis the instability domain is constructed in the 
phase space {9, </>, A, g} comprising the nonideality pa- 
rameter 9, the difference (f> of the species interaction con- 
stants, the parameter A characterizing the ratio between 
time scales of species incorporation into the crystalline 
lattice, and the species diffusion flux (in dimensionless 
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units) . The potential difference (j> > is assumed before- 
hand to be nonnegative because, otherwise, exchanging 
the species indices makes it value positive. Projection of 
this domain onto the plane {9, </>} for a fixed value of A 
enables us to divide all the points on the plane {9, </>} into 
stable and unstable ones. The latter points correspond 
to such solids for which the growth instability under con- 
sideration can arise in principle. 

It has been demonstrated that there are five charac- 
teristic regions on the plane {9, </>}, where the growth 
instability exhibits different properties. In particular, in 
the region 

{9>9 C - 9 + 0<29 cO } 

the growth instability is governed mainly by the solid 
nonideality and was analyzed in detail previously in 
Ref. [ll[. In the region 

{9>9 C ; 9 + 0>29 cO } 

for A ^ 1 the instability onset is governed by the system 
asymmetry and, as a result, only one species plays an 
active role, the diffusion flux of the other component is 
practically quasiequilibrium. However for large values of 
the diffusion flux the instability dynamics again is mainly 
affected by the solid nonideality. In the region 

{9<9 C - 9 + <j)>29 cQ } 

for A ^ 1 the instability is due to the system asym- 
metry even for large values of the diffusion flux. It can 
arise also for the ideal solid solution. In this case the 
the critical value g c of the species diffusion flux exhibits 
a rather complex behavior neat the instability boundary, 
in particular, g c remains bounded as the system comes 
close to it. It has demonstrated that the system asym- 
metry can induce, in principle, the growth instability if 
A > 1, however if A — > 1 the required value of the po- 
tential difference <fi c — ► oo (for a fixed value of 9 < 9 c q). 
The condition that the system admits an unstable per- 
turbation with finite spatial scales for large values of the 
species diffusion flux, g — > oo, gives a fairly precise ap- 
proximation of the boundary of the instability caused by 
the system asymmetry except for values of A close to its 
threshold A = 1. 

Analyzing the limits cases of the growth dynamics two 
typical regimes were singled out. One of them is the 
regime of constant diffusion flux that characterizes "slow" 



dynamics of species concentration and solid composition. 
The other referred to as the regime of constants sur- 
face concentration described the stage of "fast" dynam- 
ics. Oscillatory zoning studied in our previous paper (llj 
corresponds to the case when the region of constant sur- 
face concentration holds with respect to all the species. 
As a result the phase portrait of the system dynamics 
looks line a limit circle of relaxation oscillations on the 
phase plane {Cf/CJ,x}- At a rough approximation it 
can be constructed referring to the iV-like curve show- 
ing the quasi-stationary dependence of the ratio C{ /Cf 
on x- I n the present paper the main attention is paid 
to the case A 3> 1 where the nonlinear stage of the de- 
veloped instability is characterized by the regime of con- 
stant surface concentration with respect to one species 
and regime of constant diffusion flux with respect to the 
other species. Now the phase plane {Cf,x} gives the 
appropriate representation of the system portrait in a 
similar way, including the construction of the limit circle 
describing oscillatory zoning. 

Numerical simulation justifies these conclusions. Be- 
sides, the species distribution in the aqueous solution 
bulk found numerically demonstrates the fact that a 
rather sophisticated model of the boundary layer should 
be developed to describe oscillatory zoning adequately. 

At the next of the theory development OZ in 2D case 
will be considered with respect to two aspects. One is 
the affect of 2D species distribution itself on the pattern 
formation. The other is due to the fact that, for exam- 
ple, in the domain 2?o-i the "optimal" conditions for the 
instability development can match the solid composition 
X 7^ 0.5. In this case a special instability with respect to 
nonuniform perturbations along the crystal surface can 
arise. So we expect that the characteristic length de- 
termining spatial correlations of the OZ-pattern will be 
found in this way. Nevertheless referring to OZ obtained 
in laboratory [3, 0, 0] this length should be expected to 
exceed essentially the size of growing crystallites about 
200 /xm. 
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